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Abstract 


Persistence  in  weather  forecasting  is  used  to  describe  runs  of  sev¬ 
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sistence,  P,  are  limited  by  existence  considerations  to  —  1  <  P  <  1. 
For  P  =  0,  the  increments  making  up  the  sum  are  uncorrelated,  in¬ 
dependent  variables.  Values  of  P  near  unity  represent  a  tendency  for 
long  runs  of  similar  values.  Observed  global  average  annual  tempera¬ 
ture  records  exhibit  strong  positive  persistence  even  when  linear  trend 
is  removed.  A  hundred  year  CCM  run  (CCM-1)  shows  vanishing  per¬ 
sistence  perhaps  indicating  that  the  real  oceans  give  rise  to  runs  of 
several  years  or  decades  with  similar  climate  characteristics  while  the 
model  ocean  fixed  by  seasonal  means  does  not. 
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1  INTRODUCTION 


Questions  about  predictability  of  climate  focus  attention  on  the  phe¬ 
nomenon  of  persistence.  In  meteorology,  the  term  “persistence”  is  used  to 
describe  runs  of  several  days  with  similar  weather  characteristics.  Persistence 
can  aid  weather  forecasting;  Under  some  circumstances  today’s  weather  can 
be  used  to  predict  tomorrow’s.  In  the  case  of  high  frequency  weather,  per¬ 
sistence  is  associated  with  the  large  scale  organized  motion  within  the  at¬ 
mosphere,  the  nature  of  which  is  determined  in  large  measure  by  the  earth’s 
rotation.  For  example,  the  air  mass  comprising  a  high  pressure  system,  with 
anticyclonic  circulation,  has  similar  characteristics  over  the  area  it  covers  and 
weather  conditions  can  persist  for  days  as  the  system  moves  over  a  locality. 
Such  a  system  may  take  several  days  to  pass  over  a  given  location  in  summer, 
while  in  winter  it  moves  more  rapidly.  Therefore,  in  general,  predictive  skill 
is  greater  in  summer  than  in  winter. 

Does  climate  show  persistence?  Are  there  runs  lasting  several  years 
or  more  in  which  similar  conditions  persist  that  are  different  from  “mean” 
conditions?  Does  climate  exhibit  long  term  memory  so  that  long  term  cor¬ 
relations  in  patterns  of  climate  exist?  Are  long  runs  of  drought  recorded  in 
Biblical  times  and  in  today’s  California  examples  of  persistence  in  climate? 

If  climate  does  demonstrate  these  long  term  correlations  in  the  form 
of  runs,  then  this  persistence  may  assist  long  term  climate  prediction.  The 
physical  origin  for  such  persistence  might  be  in  the  longer  term  organized 
motion  of  the  ocean,  where  time  scales  of  a  thousand  years  may  be  expected, 
or  in  the  complexities  of  nonlinear  dynamics  of  both  oceans  and  atmospheres. 
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While  persistence  is  a  familiar  term  in  meteorology  there  is  no  quantita¬ 
tive  definition  nor  any  clear  cut  way  of  determining  the  degree  of  persistence 
in  an  observed  or  calculated  time  series.  Long  runs  of  similar  values  of  a 
parameter  describing  climate,  such  as  global  einnual  average  temperature, 
would  reveal  themselves  as  energy  in  the  low  frequency  bands  of  the  power 
spectrum.  In  fact,  long  meteorological  (and  other  geophysical)  records  typi¬ 
cally  show  abundant  energy  at  low  frequencies,  often  described  as  red  noise. 
This  behavior  is  illustrated  in  the  cumulative  spectrum  (see  Figure  1 )  of  the 
average  global  temperature  calculated  from  the  Lorenz  27-variable  model  in¬ 
tegrated  over  26,000  years  (see  Section  4).  About  90  percent  of  the  variance 
lies  in  the  lower  20  percent  of  the  frequency  domain.  Century  long  observa¬ 
tional  temperature  records  show  similar  spectra  (see  Section  5).  A  current 
GCM  with  ocean  temperatures  fixed  by  seasonal  means  does  not  display 
observed  low  frequency  behavior  nor  persistence  (see  Sections  5.2  and  5.3). 

Classical  time  series  statistics  are  weak  in  describing  low  frequency  be¬ 
havior,  although  the  two  point  correlation  function  is  useful  in  describing 
short  time  correlations.  Higher  order  correlations  are  required  to  describe 
long  term  runs,  but  these  quantities  are  very  computationally  intensive  and 
do  not  lend  themselves  to  easy  graphical  presentation  or  interpretation. 
Power  spectra  are  notorious  for  their  weakness  in  representing  behavior  near 
zero  frequency.  If  the  mean  is  removed,  the  zero  value  at  zero  frequency 
pulls  down  estimates  of  neighboring  frequencies  due  to  the  “windowing”  as¬ 
sociated  with  the  finite  length  of  the  record.  If  the  mean  is  included,  the 
delta  function  at  zero  frequency  pulls  up  nearby  estimates.  Any  linear  trend 
is  usually  removed  prior  to  taking  the  spectrum  because  the  low  frequency 
content  of  a  trend  distorts  the  estimates  of  the  low  frequency  part  of  the 
spectrum.  Since  the  number  of  “cycles”  at  periods  comparable  to  the  record 
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length  is  small,  the  statistical  reliability  of  the  low  frequency  estimates  is 
based  on  few  samples  and  is  therefore  low. 

The  low  frequency  components  in  a  time  series  can  be  accentuated  by 
filtering.  One  means  of  filtering  is  to  integrate  or,  for  equally  spaced  data, 
to  form  the  cumulative  sum.  K  /(w)  is  the  estimate  of  the  power  spectrum 
at  angular  frequency  u>  of  the  observed  time  series  x{t),  then  f[u)luP‘  is  the 
estimate  (except  at  zero  frequency)  of  the  spectrum  of  the  cumulative  sum 
X{n) 

X{n)  =  ^(0- 

t=i 

The  summing  highlights  the  low  frequency  components  of  x{t)  without  com¬ 
pletely  eliminating  high  frequency  terms,  as  would  be  the  case  for  low  pass 
filtering.  In  the  time  domain,  long  runs  of  similar  values  of  x{t)  will  be  exhib¬ 
ited  as  nearly  monotonic  increases  or  decreases  of  X(t),  while  a  series  with 
less  correlation  will  show  a  much  more  ragged  wave.  As  will  be  shown,  a 
measure  of  the  irregularity  of  the  cumulative  sum  X(n)  provides  a  measure 
of  the  persistence  within  the  calculated  or  observed  time  series  x(t). 

In  the  case  where  the  x(<)  are  normally  distributed  and  independent, 
with  vanishing  mean  and  unit  variance,  the  running  sum  X(n)  can  be  de¬ 
scribed  as  a  classical  Brownian  process.  The  increment  X(n  -i-  k)  —  X(n)  is 
a  normal  variable  with  mean  0  and  variance  k  which  is  independent  of  X(n) 
an_d  the  values  X(m)  for  m  <  n  since  x(t)  are  by  assumption  uncorrelated. 
In  a  clcissical  Brownian  process,  there  is  no  memory  of  the  past.  The  lack 
of  memory  is  a  direct  result  of  the  independent  nature  of  x(t).  Because  x(t) 
are  independent,  their  index  of  persistence,  denoted  by  P,  and  which  will 
be  defined  shortly,  will  vanish.  For  notational  purposes  we  note  that  the 
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variance  of  X(n)  scales  as  where 


i/=  1/2 

for  a  Brownian  process. 

At  the  opposite  extreme  to  a  Brownian  process  with  zero  persistence 
we  consider  a  process  in  which  the  increments  of  2:(<)  are  a  constant  random 
number  xj. 

x{t)  =  Xj 

with  E[xi]  =  0  and  E[x\\  =  1.  The  cumulative  sum  for  this  process  is  simply 

A’(n)  =  xi  n 

and  the  variance  of  vV(n)  is  given  by 

var  [X]  =  £^[xi] 

so  that  the  variance  scales  as  where  H  =  \.  For  this  process  persistence 
is  perfect,  since  the  run  length  of  identical  values  equals  the  length  of  the 
record.  If,  in  analogy  with  the  correlation  coefficient,  we  assign  P  =  1  to 
perfect  persistence  and  P  =  0  to  the  absence  of  persistence,  then  P  is  related 
to  the  scale  factor  H  by 

P  =  2H-U  -1<P<1- 

Negative  values  of  P  correspond  to  antipersistence. 

Hurst  (1951),  a  hydraulic  engineer,  first  introduced  the  use  of  cumula¬ 
tive  sums  to  analyze  long  hydrological  records  by  considering  the  storage  of 
water  in  reservoirs.  The  amount  of  water  in  a  reservoir  at  any  one  time  is 
determined  by  the  cumulative  inflow  minus  the  evaporation  and  draw  down. 
In  a  managed  reservoir,  the  inflow  fluctuates  widely  but  the  draw  down  is 
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more  or  less  constant.  Motivated  by  reservoir  considerations,  Hurst  analyzed 
a  number  of  records  using  a  rescaled  range  (R/S)  analysis  described  in  detail 
in  the  following  section.  Hurst’s  results  were  puzzling.  Hydrological  records 
did  not  scale  according  to  =  0.5  or  P  =  0  but  showed  values  of  H  that 
corresponded  to  positive  values  of  P.  For  reservoir  height,  persistence  is  un¬ 
derstandable,  since  it  depends  on  rainfall  in  prior  seasons  as  well  as  rainfall 
in  the  current  year,  and  therefore  the  system  has  memory.  Hurst  also  found 
H  >  1/2  for  long  temperature  records  where  the  origin  of  “memory”  is  not 
as  clear. 

This  paper  examines  persistence  in  climate  and  models  of  climate.  In 
section  3,  we  describe  the  theory  of  process  with  stationary  increment  and  in¬ 
troduce  the  notion  of  fractional  Brownian  motion  first  studied  by  Kolmogorov 
(1940)  and  later  by  Levy  (1953).  Such  processes  have  memory  and  positive 
persistence.  A  variety  of  methods  to  determine  P,  including  R/S  analyses, 
are  employed  in  Section  4  to  determine  the  persistence  of  long  (26,000  years) 
records  generated  by  a  low  order  climate  model.  The  results  of  applying 
these  methods  to  observational  records  are  discussed  in  Section  5. 
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2  RESCALED  RANGE  (R/S)  ANALYSIS 
OF  INDEPENDENT  RANDOM  VARI¬ 
ABLES 


1x1  order  to  illustrate  Rj S  analysis  we  first  consider  a  series  where  x{t) 
are  mutually  independent  random  variables  with  a  common  distribution  func¬ 
tion.  For  convenience  in  calculation,  the  mean  of  x(t)  is  set  to  zero  and  the 
values  are  normalized  such  that  the  variance  is  unity.  We  form  the  cumulative 
sums  of  a:(t)  in  order  to  bring  out  the  low  frequency  behavior 

X{n)  =  '£x(t). 

t=\ 


Order  statistics  are  more  robust  than  conventional  statistical  measures. 
To  apply  order  statistics  we  form  the  range  R{n)  by  setting 

M{n)  =  max  [0,X(l),A'(2),...,X(n)] 
m(n)  =  min  [0,  A'(1),A'(2),...,  A'(n)] 

and 

R{n)  =  M  (n)  —  m(n). 

For  ease  in  interpretation  it  is  useful  to  modify  this  definition.  Suppose  that 
the  total  length  of  the  record  is  N.  Instead  of  considering  the  sums  A'^(n), 
we  consider  instead  their  deviations  from  a  line  drawn  joining  the  origin  to 
the  point  [N ,X{N)).  Thus  we  replace  the  random  variables  A(n)  by 

A*(n)  =  X{n)-nX{N)IN 

and  define  the  corresponding  variables  /?*(n),  M*(n)  and  m*(n)  accordingly. 
The  random  variable  /?*(n)  will  be  called  the  adjusted  range  of  the  cumula¬ 
tive  sum. 
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In  applications,  the  adjusted  sum  is  to  be  preferred  because  it  appears 
to  have  greater  sampling  stability  (it  is  more  robust)  and  it  eliminates  a 
trend  in  case  that  £[x(/)]  ^  0.  Since  the  sums  X(N)  are  asymptotically 
normally  distributed,  the  asymptotic  distribution  of  the  range  is  independent 
of  the  distribution  of  x{t).  Summing  not  only  emphasizes  the  low  frequency 
components  of  x{t)  but  also  ensures  a  robust  quality  to  the  statistics  of  X{n), 
such  as  the  range  or  the  adjusted  range, 

A  straightforward  calculation  of  the  range  for  normally  distributed  x(t) 
leads  to 

£;[/i:(n)]  =  2(2n/jr)»/2 

and 

var  [/2(n)]  =  4n(/n2  —  2/x) 

(Feller,  1951).  The  statistics  for  the  adjusted  range  show  the  same  depen¬ 
dence  on  the  length  n  over  which  the  sum  is  taken 

£;[ir(n)]  =  (n7r/2)‘/*  ~  1.2533 
2 

var  [/2*(n)  =  (-- — —)  n  ~  0.07414  n. 

6  2, 

The  result  that  the  range  varies  as  the  square  root  of  the  length  of 
the  sum  depends  critically  on  the  assumption  that  the  variables  x[t)  are 
independent.  If  this  assumption  does  not  hold,  then  A’(n  -f-  1)  depends  on 
X{n)  ,  and  since  X{n)  is  determined  by  the  previous  n  values  of  x{t),  A'^(n-l-l) 
depends  on  x(<),  t  =  l...n.  The  cumulative  sums  provide  a  measure  of  past 
correlation  if  the  dependence  of  the  range,  or  the  adjusted  range,  does  not 
vary  as  the  square  root  of  the  summation  interval. 

The  quantity  Rf S  is  formed  by  dividing  the  adjusted  range  of  A(n) 
by  5(n),  which  is  the  standard  deviation  of  z(n).  For  stationary  Gaussian 
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processes,  division  by  S{n)  adds  little  to  the  analysis;  x{t)  can  always  be 
normalized  to  have  unit  variance.  However,  for  nonstationary  processes  with 
long  range  dependence,  division  by  S{n)  leads  to  a  more  stable  statistic. 

Based  on  the  above  considerations,  we  can  apply  R/S  analysis  to  test 
for  correlation  in  random  number  generators  commonly  employed  in  com¬ 
puter  packages.  The  random  numbers  were  generated  according  to  the  algo¬ 
rithm  proposed  by  Park  and  Miller  (1988),  which  is  under  consideration  for 
adoption  as  a  IEEE  standard.  Figure  2  shows  a  sequence  of  8192  normally 
distributed  variates  supposed  to  be  uncorrelated.  The  observed  distribu¬ 
tion  (see  Figure  3)  closely  approximates  the  normal  distribution.  Figure  4 
illustrates  the  cumulative  sum  of  the  variables  in  Figure  2  with  the  ragged 
curve  typical  of  a  Brownian  process.  The  R/S  statistics  are  obtained  by  first 
calculating  the  adjusted  range  i2“  for  the  total  record  of  length  N  =  8192  and 
normalizing  by  the  standard  deviation  of  the  record  of  length  N.  The  pro¬ 
cess  is  repeated  for  chunks  of  record  of  length  N/2,N/4,...N/102A,  and  the 
resulting  value  of  the  scaled  adjusted  range,  R/S,  is  obtained  by  averaging 
over  the  2,4,  ...1024  non-overlapping  intervals.  The  smallest  interval  contains 
eight  points.  The  concept  of  range  begins  to  break  down  at  small  interval 
and  chunks  smaller  than  8  are  not  used.  A  plot  of  the  log  {R/S)  against 
log  (n),  where  n  ’s  the  length  of  the  interval  over  which  R/S  is  estimated, 
determines  the  scaling  exponent  H  in  the  relation 

R/S  n^. 

For  the  series  shown  in  Figure  2  the  scaling  slope  is  0.51  (see  Figure  5).  For 
n  =  8192,  the  difference  in  adjusted  range  for  slopes  0.5  and  0.51  is  11.8, 
while  the  standard  deviation  of  the  estimate  of  adjusted  range  is  24. o9,  so 
that  the  value  0.51  is  not  statistically  distinct  from  0.5.  For  comparison, 
a  slope  of  0.6  would  produce  an  adjusted  range  more  than  eight  standard 
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Figure  2.  Time  sequence  of  8192  normaUy  distributed  random  variables  with  zero  mean  and 
unit  variance. 
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deviations  from  the  expected  value  for  a  random  independent  variable. 

We  note  that  while  the  range  for  the  cumulative  sum  of  the  normally 
distributed  independent  variables  scales  with  length  of  observation  as 
the  range  of  the  elements  in  the  sum  x{t)  increases  feir  more  slowly,  essentially 
as  (log  (Leadbetter  et  al.,  1980;  Abarbanel  et  al.,  1991). 

A  trend  introduces  a  long  term  persistence  into  a  record,  as  would  a  pure 
harmonic,  A  sin(27r f  +  <i>).  We  take  the  frequency  /  and  amplitude  A  as  fixed 
but  the  phase  4>  is  taken  as  random  having  a  uniform  distribution.  For  x{t) 
equal  to  a  sine  wave,  the  range  of  the  cumulative  sum  X{t)  in  the  limit  of 
large  n  is  A/irf  while  the  limit  of  the  standard  deviation  is  A/2,  so  that  the 
limit  of  R/ S  is  2/7r/.  The  statistic  R/S  for  a  pure  sine  wave  does  not  vary 
according  to  so  it  shows  persistence  (actually  perfect  antipersistence). 

The  situation  differs  in  the  limit  of  large  n  when  white  Gaussian  noise 
of  zero  mean  and  unit  variance  is  added  to  the  process 

i(t)  =  Asin(27r/<  4>)-\-  g{t). 

The  cumulative  sum  A(n)  then  satisfies  the  inequality 

G'(n)  -  A/2irf  <  XCn)  <  G(n)  +  A/27r/ 
where  G(n)  is  the  cumulative  sum  of  the  noise  g(t). 

In  the  limit  of  large  n,  the  ranges  of  X(n)  and  G(n)  are  identical  since 
the  range  of  X(n)  varies  as  The  variance  of  the  sine  wave  plus  noise  in 
the  limit  of  large  n  is 

ElS(n)J  =  1  +  A/2 
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and  therefore  the  statistic  RjS  h«is  an  expected  value  for  large  n  of 

The  value  of  RJS,  as  well  as  the  time  required  to  reach  the  limit,  is 
determined  by  Al2ir f.  Very  low  frequency  sine  waves  will  lead  to  persistence 
in  the  presence  of  uncorrelated  noise.  Lines  in  the  spectrum  at  low  frequencies 
will  affect  the  estimate  of  R/S,  as  will  a  linear  trend.  High  frequencies  will 
not  show  up  in  the  R/S  statistic. 

The  existence  of  long  term  correlations  can  have  important  consequences. 
For  example,  the  physical  model  for  Brownian  motion  is  one  of  a  small  par¬ 
ticle  being  bombarded  by  molecules  in  thermal  equilibrium.  The  depen¬ 
dence  flows  from  the  assumption  that  the  molecules  in  thermal  equilibrium 
have  uncorrelated  motions.  With  vanishing  correlation  the  ordinary  diffusion 
equation  applies.  Statistical  independence  at  large  time  and/or  space  scales 
is  an  essential  ingredient  of  the  concept  of  thermal  equilibrium;  the  exis¬ 
tence  of  longer  term  correlations  would  invalidate  the  assumption  of  thermal 
equilibrium. 
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3  FRACTIONAL  BROWNIAN  MOTION  AND 
PERSISTENCE 


For  the  moment  we  return  to  the  cumulative  sum  of  independent,  sta- 

n 

tionary  increments,  and  form  the  structure  function  for  X{k),  X{k)  = 

t=i 

defined  by 

Z)(A:)  =  £[(X(A:)-X(0))^] 

where  by  convention  we  can  set  A’(O)  =  0.  The  structure  function  for 
X{n  +  m)  is 

D{n  +  m)  =  f:[(X(n  +  m)-X(0)f]  =  F:[(X(n  +  m)-X(m))'] 


+  £:[(X(m)-X(0))% 


But  we  also  have 


E[{X{n  +  m)  -  X(m))^]  =  D{n) 


if  the  increments  are  stationary  and  independent  and 


D{n  +  m)  =  D{n)  +  D{m)  for  n  >  0,  m  >  0  . 


The  only  possible  functional  form  for  the  structure  function  is 


D{n)  =  Cn, 

identical  to  the  functional  form  for  the  variance  of  the  range  described  in  the 
previous  section. 


The  spectral  density  F(w)  associated  with  X{t)  is  related  to  the  contin¬ 
uous  structure  function  D{t)  by 


D{t)  =  4  /  (1  —  cos  u  t)  F(uj)  du>  . 

Jo 
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The  spectrum  of  the  independent,  random  increments  is  a  constant 


therefore  the  spectrum  of  the  cumulative  sum  of  independent,  random  incre¬ 


ments  is 


Fiu;) 


C 

2nu)^' 


Kolmogorov  (1940)  generalized  the  concept  of  Brownian  processes  to 
consider  structure  functions  of  the  form 

D{t)  = 

with  0  <  H  <  1.  The  limits  on  H  come  from  examining  the  spectral  repre¬ 
sentation  of  D{t).  Since  from  dimensional  considerations 


and  since  1  —  cos  u;r  varies  as  near  zero  frequency,  the  spectral  represen¬ 
tation  for  D{t)  exists  if  2H  +  1  <  3,  which  fixes  the  upper  limit  for  H.  At 
high  frequencies,  an  ultraviolet  catastrophe  can  only  be  avoided  if  //  >  0. 
These  limits  on  H  allow  the  persistence  index  P  =  2H  —  1  to  run  from  -1  to 
1.  Since  the  spectrum  for  X  can  be  written  as 


c 


the  proportionality  constant  Ci  is  related  to  C  in  the  relation  for  the  structure 
function  by 

C  =  ^ 

'  4  /o  (1  —  cos  z)  x“2^~*dz 


ri2H  +  l)sm{2HT/2)C 
2v 
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The  power  law  dependence  of  the  structure  function  implies  that  D  and 
X  are  self  similar.  Thus  a  change  in  scale  in  the  time  variable 

n  =  hn 


leads  to  a  scale  change  in  X 


X{n)  =  aih)X{hn) 


For  classical  Brownian  motion  we  have 


a  = 


or  more  generally 


a  =  h 


-H 


and  correspondingly  the  variance  scales  as 


var[A’(n)]  =  var  [A'’(/in)] 


The  self  similar  character  of  the  cumulative  sums  implies  that  the  “box 
dimension”  of  the  curve  is  not  unity.  The  box  dimension  is  determined  by 
the  relation  between  the  number  of  boxes  required  to  cover  the  curve  and 
the  dimension  of  the  individual  boxes.  If  the  length  of  the  record  is  N ,  then 
NIhn  segments  of  length  hn  are  needed  to  cover  the  time  axis.  In  each 
segment  the  range  of  the  record  is  R{hn)  =  R{n)  and  we  need  a  stack 
of  R{n) / ha  boxes  of  height  ha  to  cover  the  range.  The  total  number  of 
boxes  is  then 


R{n)  T  iH-2  i-Da 

B{h;a,n)  =  — 7 - —  ~  ~ ® 

ha  hn 

where  Db  is  defined  as  the  box  dimension.  In  terms  of  the  parameter  H  the 
box  dimension  is 


Db  =  2-H. 
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The  determination  of  the  box  dimension  for  the  sequence  of  random 
variable  shown  in  Figure  2  is  illustrated  in  Figure  6.  Here  the  box  dimension 
Db  =  1-48,  which  yi'^lds  an  H  value  of  0.52;  this  can  be  compared  to  the 
value  of  0.51  obtained  by  an  RjS  analysis. 

Series  with  values  of  the  parameter  H  >  \ /2,  or  P  >  0,  show  persistence; 
that  is,  the  process  generating  the  time  series  has  memory.  Levy  (1953) 
demonstrated  this  by  showing  that  for  X{t)  generated  as  a  Brownian-like 
process,  X{t)  can  be  represented  as 

where  dVF(s)  is  the  increment  of  an  ordinary  Brownian  process  {H  =  1/2). 
This  representation  clearly  demonstrates  the  long  term  memory  associated 
with  a  fractional  Brownian  process. 


Persistence  of  a  fractional  Brownian  process  can  also  be  shown  by  ex¬ 
amining  the  correlation  properties.  In  particular,  the  past  increments  are 
correlated  with  future  increments.  Given  the  increment  A'(O)  —  X(— n)  from 
time  n  to  time  0,  the  probability  of  having  an  increment  X(n)  —  A^(0)  aver¬ 
aged  over  the  distribution  of  past  increments  is  proportional  to 

£[(A:(0)  -  A-(-n))(A'(n)  -  A-(O))]. 


For  convenience  we  choose  the  origin  so  that  X(0)  =  0.  The  correlation 
function  of  future  increments  with  past  increments  is  then 

£:[_X(-n)X(n)] 


C(n)  = 


E[X{n)^] 


where  we  have  normalized  with  the  variance.  The  above  scaling  relations 
£[(A:(n)  -  X{kY]  ~  |n  -  lead  to 


C(n)  =  2^'^'*  -1  =  2^-1. 
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Figure  6.  Determination  of  box  dimension  of  the  cumulative  sum  curve  shown  in 
Figure  4  for  a  zero  mean  unit  variance  random  variable  having  a  normal 
distribution  Slope  =  -1.478S.  which  corresponds  to  H  =  0.521 1. 
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The  correlation  of  future  increments  with  past  increments  vanishes  for  strict 
Brownian  motion.  However,  for  positive  persistence,  the  correlation  is  posi¬ 
tive  and  independent  of  the  time.  This  implies  an  increcising  trend  in  the  past 
will  produce  an  increasing  trend  in  the  future.  In  fact,  in  time  limited  sys¬ 
tems,  fractional  Brownian  motion  is  an  approximate  model  over  some  range 
of  time  scales  but  not  necessarily  for  the  entire  record. 

Fractional  Brownian  motion  is  illustrated  in  Figures  7  and  8.  In  Figure 
7,  the  curve  corresponds  to  a  persistence  of  0.8.  The  long  runs  of  x{t)  with 
the  same  sign  show  up  clearly  as  nearly  straight  climbs  or  descents  in  the 
cumulative  sum  A'^(n).  Antipersistence  is  shown  in  Figure  8,  giving  rise  to  a 
highly  ragged  curve. 
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Figure  7.  Fractional  Brownian  motion,  P  =  0.8 
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4  ANALYSIS  OF  LONG-TIME  INTEGRA¬ 
TION  OF  LORENZ  27- VARIABLE  MODEL 


4.1  Introduction 


We  use  values  generated  by  a  Lorenz  27-variable  model  of  the  atmo¬ 
sphere  to  illustrate  the  analysis  described  in  Sections  2  and  3.  Lorenz  (1984) 
formulated  a  “low  order”  model  of  atmospheric  change  with  a  goal  of  explor¬ 
ing  a  “moist”  atmosphere.  A  moist  atmosphere  contains  bulk  liquid  water 
and  water  vapor  that  can  form  clouds.  A  shallow  ocean  provides  water  for 
evaporation  as  well  as  a  heat  sink  or  source.  Clouds  form  or  dissipate  in 
response  to  changes  in  relative  humidity  with  consequent  changes  in  albedo 
and  the  radiative  balance. 

In  the  model,  the  albedo  is  proportional  to  the  cloud  cover;  since  freezing 
is  excluded  from  the  model  there  are  no  changes  in  albedo  due  to  varying 
snow  and  ice  cover.  The  cloud  cover  is  parameterized  in  terms  of  the  relative 
humidity.  Locally,  the  relative  humidity  drops  with  higher  temperature  and 
the  albedo  decreases,  leading  to  more  solar  radiation  to  heat  the  earth.  The 
cloud-albedo  feedback  allows  for  significant  swings  in  local  and  global  average 
temperature. 

The  radiative  and  thermodynamic  processes  that  characterize  a  moist 
atmosphere  introduce  complicated  nonlinear  terms  into  the  governing  equa¬ 
tions.  In  order  to  limit  the  computational  complexity,  Lorenz  (1984)  intro¬ 
duced  a  number  of  simplifying  approximations.  Lorenz  found  that  for  short 
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term  integrations  the  model  produced  results  that  qualitatively  were  in  rea¬ 
sonable  agreement  with  observed  atmospheric  behavior  and  concluded  that 
the  model  was  suitable  for  production  runs.  Abarbanel  et  al.  (1990)  and 
Abarbanel  et  al.  (1991b)  investigated  the  climate  of  the  Lorenz  model  by 
integrations  over  several  thousand  years.  We  have  continued  these  studies 
using  a  25,521  year  integration  with  a  1.25  hour  time  slip  and  averaging  over 
the  year  and  over  the  surface  of  the  model  to  obtain  “global”  annual  average 
temperature.  The  model  does  not  contain  either  a  seasonal  or  diurnal  cycle. 

An  important  point  to  note  is  that  the  model  contains  only  one  long 
time  constant:  the  thermal  time  constant  for  warming  the  shallow  ocean, 
about  ten  years.  The  ocean  is  taken  to  have  ten  times  the  heat  capacity  of 
the  atmosphere.  Horizontal  heat  flow  in  the  ocean  is  not  allowed.  Thus  there 
is  a  permanent  equator  to  pole  temperature  gradient  that  undergoes  small 
fluctuations  in  response  to  the  poleward  transfer  of  heat  by  the  atmosphere. 
The  pole  to  equator  thermal  gradient  forms  a  large  scale  coherent  structure 
with  a  time  scale  equal  to  the  length  of  the  record. 

4.2  Statistics  of  the  Lorenz-27  Model 


The  results  of  a  25,521  year  integration  of  Lorenz-27  for  global  annual 
average  temperature  are  d’^'played  in  Figure  9.  The  temperatures  have  been 
scaled  to  unit  variance  with  a  standard  deviation  of  0.364°K  (see  Table  1). 
The  time  history  of  the  temperature  resembles  a  white  noise  process  and  any 
long  term  trends  or  persistence  are  m^lsked  by  the  irregular  short  term  behav¬ 
ior.  The  histogram  for  the  temperature  values  is  shown  in  Figure  10.  The 
temperature  values  closely  approximate  a  normal  distribution.  The  skew- 
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Time  in  Years  (sd  =  0.36  "K) 


Figure  9.  Variation  in  global  annual  average  temperature  for  a  25,521  year  integration  of  the 
Lorenz  27-variable  model  The  values  are  reduced  to  zero  mean  and  unit  variance. 
The  standard  deviation  of  the  raw  temperature  record  is  0.364 °K. 
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Table  1 


Standard  Deviation  of  Global  Annual  Average 
Temperature  of  Long  (25,521  years)  Integration 
of  Lorenz  27- Variable  Model 


Interval 

(Years) 

Standard  Deviation 
(Degrees  K) 

1  -  25,521 

0.3.64 

50001  -  15  000 

0.374 

15001  -  25000 

0.359 

40003  -  4302 

0.270 

23401  -  23700 

0.321 

12374  -  12673 

0.342 

14001  -  14140 

0.341 

ness  coefficient  is  0.18,  which  is  not  significant.  The  distribution  also  closely 
approximates  a  normal  distribution  in  the  tails.  The  range  of  the  tempera¬ 
ture  in  units  of  standard  deviation  is  shown  in  Figure  11  together  with  the 
expected  value  for  a  variable  drawn  from  a  normal  population.  The  auto¬ 
correlation  function  for  the  temperature  record  is  shown  in  Figures  12  and 
13.  The  correlation  shows  an  exponential  decay  with  a  first  zero  at  about 
65  years.  At  longer  lags  the  correlation  function  shows  structure,  though  the 
paucity  of  samples  for  long  lags  raises  question  about  statistical  significance. 

The  power  spectrum  of  Lorenz-27  shows  abundant  energy  at  low  fre¬ 
quencies  (see  Figures  14  and  15).  The  spectrum  shows  a  number  of  peaks 
at  periods  greater  than  100  years.  The  cumulative  spectrum  (Figure  1) 
establishes  that  90  percent  of  the  variance  of  the  record  displayed  in  Fig¬ 
ure  9  is  contained  at  frequencies  less  than  0.1  cycles  per  year  (cpy).  For 
periods  less  than  16  years,  the  spectrum  drops  off  as  /“2.2  Figure  16). 
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Figure  11. 
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Variation  in  range  of  the  temperature  values  shown  in  Figure  9.  The  dashed  curve 
represents  the  expected  value  of  the  range  for  a  normally  distributed  random 
variable.  (Abarbanel  et  al.  1991a) 
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Autocorrelation 
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Figtiic  12.  The  autoconelation  function  for  the  time  series  shown  in  Figure  9.  The  correlation 
function  obtained  by  taking  the  Fourier  transform  of  the  power  spectrum  is  shown 
in  Figure  14. 
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Autocorrelation 


Power 


Figure  15.  The  low  frequency  portion  of  the  power  spectrum  shown  in  Figure  14 
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Log  Frequency  (cpy) 


Figure  16.  The  high  frequency  portion  of  the  spectrum  given  in  Figure  14  in  a  log-log 

representation.  The  best  least  squares  line  fit  to  the  spectrum  is  f  ^  ^  suggesting  that 
for  time  scales  of  16  years  or  less,  the  record  given  in  Figure  9  can  be  modeled  as 
a  Brownian  motion  with  independent  increments. 


Thus  at  high  frequencies  the  record  resembles  a  classical  Brownian  walk 
derived  from  white  noise.  The  behavior  deviates  greatly  from  this  model  at 
low  frequencies,  which  contain  the  bulk  of  the  variance. 


4.3  R/S  Analyses  of  Lorenz-27  Model 


An  R/S  analysis  of  the  record  shown  in  Figure  9  yields  a  persistence 
index  P  —  0.34  for  intervals  greater  than  about  20  years  (see  Figure  17).  An 
alternative  interpretation  of  the  analysis  is  shown  in  Figure  18.  In  this  in¬ 
terpretation  there  are  three  scales  for  the  variation  in  temperature.  At  high 
frequencies,  periods  smaller  than  20  years,  the  variations  can  be  modeled 
as  classical  Brownian  motion  with  independent  increments.  For  periods  be¬ 
tween  20  years  and  400  years,  the  record  has  high  persistence,  P  =  0.64.  At 
still  longer  intervals  (400  years  to  25,000  years),  the  record  shows  weak  per¬ 
sistence,  P  =  0.07.  Whatever  correlation  exists  at  intermediate  time  scales 
is  lost  at  long  (>  400  years)  time  scales. 

In  terms  of  the  spectrum  of  the  record,  in  the  frequency  interval  cor¬ 
responding  to  time  intervals  between  400  and  20  years,  the  power  spectrum 
should  drop  off  as  .  The  observed  drop  off  yields  an  estimate  for  P  of 
0.71  (see  Figure  19)  as  opposed  to  the  estimate  of  0.64  derived  from  the  R/S 
analyses.  The  spectral  estimate  of  P  depends  sensitively  on  the  highest  fre¬ 
quency  used  in  obtaining  the  slope,  since  at  high  frequencies  the  slope  of  the 
spectrum  is  -2.2  (see  Figure  16).  Persistence  is  closely  tied  to  the  character 
of  the  spectrum  at  low  frequencies.  The  higher  the  fraction  of  the  variance 
that  is  carried  by  the  low  frequencies,  the  more  likely  that  a  record  will  show 
persistence.  In  order  to  distinguish  persistence  due  to  low  frequency  lines  in 
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L0g  (R/S) 


Log  Interval 


Figure  17.  Results  of  an  R/S  analysis  of  the  record  given  in  Figure  9.  The  logarithm  of  the 

scaled  range  is  plotted  against  the  logarithm  of  the  interval  over  which  the  rescaled 
range  was  estimated. 
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Log  (R/S) 


Figure  18.  Alternative  intepretation  ot  the  R/S  analyses  of  the  record  shown  in  Figure  9. 
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Log  of  Power 


the  spectrum  from  persistence  due  to  a  continuous  low  frequency  spectrum 
a  much  longer  record  would  be  required  (see  Section  2). 

The  Lorenz-27  climate  model  shows  three  time  intervals  with  differing 
behavior.  At  high  frequencies,  the  model  behaves  as  an  oscillator  excited 
by  white  noise.  At  an  intermittent  range  of  frequencies,  the  variations  in 
global  annual  average  temperature  show  partial  coherence  as  measured  by 
the  persistence  index.  Over  these  time  scales  one  may  be  able  to  make 
predictions  about  average  behavior  provided  the  averages  are  longer  than  20 
years.  At  still  longer  periods  (greater  than  400  years)  partial  coherence  is 
h-ffp  y  lost  and  attempts  at  prediction  might  be  expected  to  fail. 


4.4  Analysis  of  Persistence  in  Short  Records 


The  estimates  of  persistence  discussed  in  Section  4.3  were  based  on  the 
analysis  of  the  entire  record.  In  climate  studies  good  global  annual  average 
temperature  are  very  much  limited  in  length,  to  about  110  to  140  years, 
though  there  is  a  simple  station  record  covering  some  300  years  (Abarbanel 
et  al.  1991b).  Table  2  lists  R/S  estimates  of  persistence  for  shorter  intervals 
of  the  record.  The  higher  values  of  persistence  shown  for  randomly  selected 
records  of  length  140  and  300  years  are  due  to  persistence  for  time  intervals 
of  20  to  400  years.  The  shorter  records  do  not  sample  intervals  greater  than 
400  years  and  thus  do  not  pick  up  the  lack  of  persistence  for  intervals  longer 
than  400  years. 
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Table  2 


Persistence  in  Loren2-27  Model  as  Estimated 
from  Intervals  of  Various  Lengths 


Interval 

Length  (Years) 

Persistence 

All  Intervals  20-400  Years 

1  -  25521 

25521 

0.34  0.64 

5001  -  15000 

10000 

0.42  0.74 

12374  -  12673 

300 

0.56 

4003  -  4302 

300 

0.61 

23401  -  23700 

300 

0.50 

14001  -  14140 

140 

0.48 

4.5  Alternative  Methods  of  Estimating  Persistence 


The  Rj S  method  of  estimating  persistence  appears  robust  since  it  de¬ 
pends  on  use  of  order  statistics  (David,  1981;  Huber,  1972).  Two  alternative 
methods  to  R/S  statistics  are  the  analysis  of  the  variation  of  the  variance 
with  length  of  the  interval  and  the  estimation  of  the  box  dimension  (see 
Section  3).  In  the  variance  method,  the  variance  of  the  cumulative  sum  is 
calculated  for  the  whole  record  and  then  for  successive  smaller  non  overlap¬ 
ping  intervals  obtained  by  repeated  halving  of  the  length  of  the  intervals. 
The  dependence  of  the  variance  on  the  length  of  the  interval  then  determines 
P  -f  1,  where  P  is  the  persistence  index. 

Table  3  provides  comparison  of  various  estimate  of  the  index  of  per¬ 
sistence.  The  R/S  estimate  and  the  estimate  of  the  box  dimension  of  the 
cumulative  sum  of  the  record  are  close.  The  estimate  obtained  by  examining 
the  dependence  of  the  variance  on  interval  length  consistently  gives  a  lower 
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Table  3 


Comparison  of  Various  Methods 
of  Estimating  Persistence 


Record 

R/S 

Variance 

Box  Dimension 

Lorenz  1  -  25521 

0.34 

0.22 

0.31 

Lorenz  5001-15000 

0.42 

0.30 

0.44 

Random  1-8192 

0.02 

-0.02 

0.01 

value  of  P  than  that  provided  by  R/S  analysis  or  the  box  dimensions. 


4.6  Comments  on  Estimating  Persistence 


The  principal  requirement  for  estimating  persistence  is  a  long  record. 
As  discussed  above,  the  value  of  P  depends  on  the  length  of  the  record, 
particularly  if  persistence  becomes  small  for  intervals  longer  than  some  in¬ 
terval.  The  value  of  P  obtained  for  short  records  should  be  viewed  with 
caution.  The  R/S  estimate  exhibits  greater  statistical  stability  than  does 
the  variance  estimate. 
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5  PERSISTENCE  IN  GLOBAL  ANNUAL 
AVERAGE  TEMPERATURE 


5.1  Observed  Temperature  Records 


The  record  of  global  annual  average  temperature  prepared  by  Jones  et 
al.  (1986  a,b)  and  brought  up  to  date  by  Jones  (1988)  and  Jones  and  Wigley 
(1990)  is  shown  in  Figure  20.  The  record  combines  observations  taken  both 
at  land  stations  and  ships  at  sea.  The  records  are  fraught  with  problems 
related  to  the  homogeneity  of  land  and  marine  data,  but  represent  the  best 
information  available.  The  values  earlier  than  1890  are  undoubtedly  less 
reliable  than  those  in  subsequent  years. 

The  temperature  time  series  in  Figure  20  shows  a  trend  with  warming 
apparent,  particularly  in  the  decade  of  the  1980s.  A  leaist  squares  fit  of  a 
linear  trend  provides  a  slope  of  0.028°K/year.  If  the  linear  trend  is  removed 
the  remaining  time  series  hats  a  histogram  that  approximates  a  normal  dis¬ 
tribution  (see  Figure  21).  The  standard  deviation  for  the  series  including  the 
trend  is  0.17°K;  with  the  trend  removed,  the  standard  deviation  is  0.13°K. 
The  variations  of  the  range  of  temperature  with  time  provide  a  mea-sure  of 
the  nature  of  the  tails  of  the  distribution.  As  illustrated  in  Figure  22,  the 
observed  range  approximates  the  expected  value  of  the  range  of  a  normally 
distributed  random  variable. 

Figure  23  shows  the  power  spectrum  of  the  global  annual  average  tem¬ 
perature.  The  power  is  concentrated  at  the  lower  frequencies  with  60  percent 
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Temperature  (C) 


Temperature  (sd) 
sd  =  0,13  'K 


Figure  21.  Histogram  for  the  values  of  temperature  shown  in  Figure  20  with  a  linear  trend 
removed  The  standard  deviation  for  the  record  is  O  IB^K. 
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Range 


Time  (years) 


Figure  22.  Variation  of  the  range  of  temperature  time  series  of  Figure  20  from  which  a  linear 
trend  has  been  removed.  The  dashed  curve  represents  the  expected  value  of  the 
range  for  a  normally  distributed  random  variate  with  zero  mean  and  unit  variance. 
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Power 


Power  Spectrum 


Figure  23.  Power  spectrum  of  record  shown  in  Figure  20  after  removal  of  a  linear  trend. 
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Table  4 


Estimates  of  Persistence  Obtained  for 
Global  Annual  Average  Temperature  Records 


Length 

Standard  Deviation 

Record 

(Years) 

(“K) 

Index  of  Persistence 

Hansen  and 

108 

0.22 

0.88 

Lebedeff  (1988) 
sarp#*  with  bnoar 

108 

0.14 

0.71 

trend  removed 

Jones  et  al.  (1986b) 

134 

0.17 

0.78 

same  with  linear 

134 

0.13 

0.56 

trend  removed 

Manley  (1974) 

318 

0.68 

0.46 

same  with  linear 

318 

0.61 

0.38 

trend  removed 

of  the  variance  frequencies  lower  than  0.1  cycles  per  year  (see  Figure  24).  The 
fraction  of  low  frequency  variance  is  less  than  that  displayed  in  the  Lorenz 
model  (90  percent),  but  the  observed  record  certainly  contains  uncorrelated 
measurement  error  that  may  well  reduce  the  fraction  of  energy  due  to  low 
frequency  natural  variability. 

The  persistence  index  obtained  from  an  RjS  analysis  of  the  record 
shown  in  Figure  20  is  0.78.  If  the  trend  is  removed  the  persistence  drops 
to  0  56  (see  Figure  25).  The  persistence  is  comparable  to  that  shown  in  the 
Lorenz  model  for  short  intervals  (see  Table  2).  Table  4  gives  the  value  of 
the  persistence  as  given  by  RjS  analysis  of  two  other  long  term  temperature 
records.  Hansen  and  Lebedeff  (1988)  reduced  temperature  observations  for 
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Variance 


Frequency  (cpy) 


Figure  24.  Cumulative  spectrum  of  the  time  series  given  in  Figure  20 


Log  {R  S) 


land  stations  only  to  obtain  a  global  annual  mean  surface  temperature.  The 
Manley  (1974)  record  is  pieced  together  from  observations  taken  at  various 
stations  in  Central  England  and  is  the  longest  available  instrumented  temper¬ 
ature  record.  All  the  records  show  evidence  for  significant  persistence.  The 
lower  persistence  fo:  the  Manley  record  is  probably  due  to  the  much  higher 
noise  level  associated  with  a  single  station  and  with  primitive  observations 
during  the  early  part  of  the  record. 

5.2  Persistence  in  a  Global  Circulation  Model  (GCM) 


Through  the  efforts  of  Robert  Chervin  we  had  available  a  101  year  run 
of  the  Community  Climate  Model-1  (CCM-1),  a  widely  used  general  circula¬ 
tion  model  (GCM).  This  model  contains  far  greater  complexity  than  does  the 
Lorenz  27-variable  model.  Clouds,  topography  and  snow  and  ice  feedback  are 
all  included.  The  oceans  exchange  heat  and  moisture  with  the  atmosphere 
but  the  temperature  ol  the  ocean  varies  according  to  the  seasonal  climato¬ 
logical  means.  There  are  no  ocean  currents  and,  like  the  Lorenz  model,  there 
is  no  horizontal  heat  transport  within  the  ocean.  The  resulting  global  annual 
mean  for  *he  model  is  shown  in  Figure  26.  The  standard  deviation  of  0.043°K 
is  small  compared  wnth  the  observed  records  and  with  the  Lorenz  model  (see 
Table  4). 

The  gross  statistics  of  CCM-1  are  indicated  in  Figures  27  and  28. 
About  the  mean,  CCM-1  approximates  a  normally  distributed  variate  but 
the  observed  range  is  significantly  greater  than  that  expected  for  a  normally 
distributed  variate.  The  power  spectrum  for  the  model  results  differs  dra¬ 
matically  from  the  observed  records  and  the  spectra  of  the  Lorenz  model 
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Temperature  (sd) 


20  40  60  80  100 

Time  (years) 
sd  =  0  426=K 


Figure  26.  Global  annual  average  temperature  for  a  101  year  record  produced  by  Community 
Climate  Model-1  (Chervin,  1991  personal  communication) 
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Figure  27.  Histogram  for  the  temperature  values  shown  in  Figure  26. 
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Range 


Time  (years) 


Figure  28.  Variations  of  range  for  the  temperature  series  given  in  Figure  26.  The  dashed  curve 
provides  the  expected  value  of  the  range  for  a  normally  distributed  random 
variable. 
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(see  Figures  29  and  30).  The  CCM-1  results  do  not  show  the  red  spec¬ 
trum  characteristic  of  observed  temperature  records  and  exhibited  by  the 
Lorenz  model.  In  part  the  difference  arises  from  the  treatment  of  the  oceans. 
The  energetics  of  the  ocean  determines  the  “slow”  physics.  The  prescrip¬ 
tion  of  the  boundary  condition  over  the  ocean  as  the  climatological  mean 
ensures  that  the  model  reproduces  seasonal  fluctuations  approximately  cor¬ 
rectly.  The  CCM-1  representation  does  not  exhibit  longer  term  and  higher 
amplitude  fluctuations  seen  in  nature.  Since  the  ultimate  goal  of  models  is  to 
“predict”  long  term  changes  of  climate  flowing  from  anthropogenic  composi¬ 
tional  changes,  the  lack  of  faithful  representation  of  slow  changes  is  a  major 
deficiency. 

The  101-year  run  of  CCM-1  shows  a  persistence  index  between  0.1  and 
0.2  depending  on  the  method  used  to  estimate  persistence.  These  low  values 
and  the  shortness  of  the  record  raise  questions  as  to  whether  the  estimated 
persistence  is  statistically  significant. 


5.3  Concluding  Observations 


The  analysis  presented  above  illustrates  the  usefulness  of  having  a  single 
number  P,  the  index  of  persistence,  characterize  the  long  term  behavior  of 
a  time  series.  The  existence  of  a  large  positive  value  of  P  over  a  range  of 
time  scales  indicates  that  predictability  is  possible  in  the  sense  that  runs 
or  trends  are  likely  to  be  preserved  over  these  time  intervals.  Long  term 
trends  such  as  that  due  to  changing  atmospheric  composition  give  rise  to 
persistence.  Sharp  lines  in  the  low  frequency  part  of  the  spectrum  also  result 
in  persistence.  On  short  time  scales,  persistence  is  due  to  the  existence  of 
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ince 


Cumulative  Spectrum 


Frequency  (cpy) 


Figure  30.  Cumulative  spectrum  for  global  annual  average  temperature  obtained  for  Lorenz- 
27  model,  the  observed  temperature  (see  Figure  21)  and  the  GCM  results  given  in 
Figure  26 


large  coherent  structures  in  the  atmosphere.  In  the  case  of  climate,  large  scale 
coherent  structures  in  the  ocean  may  give  rise  to  persistence  over  decades  or 
longer. 

The  analysis  of  observed  global  annual  average  temperature  records 
shows  that  these  records  have  a  strong  persistence.  While  the  records  are 
short,  only  about  a  century,  the  persistence  is  highly  statistically  significant. 
The  analysis  of  one  GCM  record  of  101  years  does  not  show  persistence  or 
a  strong  red  spectrum  and  thus  differs  significantly  from  the  observed  global 
annual  average  temperature  records.  The  lack  of  energ>'  at  low  frequencies  is 
most  probably  due  to  a  failure  to  deal  with  the  oceans  in  an  adequate  way. 
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